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HIGHER ANALOGUES OF THE DISCRETE-TIME TODA 
EQUATION AND THE QUOTIENT-DIFFERENCE ALGORITHM 
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Abstract. The discrete-time Toda equation arises as a universal equation 
^_^ ' for the relevant Hankel determinants associated with one-variable orthogonal 

Cn polynomials through the mechanism of adjacency, which amounts to the inclu- 

sion of shifted weight functions in the orthogonality condition. In this paper wc 
fl\ • extend this mechanism to a new class of two-variable orthogonal polynomials 

["f , ' where the variables are related via an elliptic curve. This leads to a 'Higher or- 

der Analogue of the Discrete-time Toda' (HADT) equation for the associated 
Hankel determinants, together with its Lax pair, which is derived from the 
Ch ' relevant recurrence relations for the orthogonal polynomials. In a similar way 

as the quotient-difference (QD) algorithm is related to the discrete-time Toda 

equation, a novel quotient-quotient-difference (QQD) scheme is presented for 

r^ , the HADT equation. We show that for both the HADT equation and the QQD 

• ■ scheme, there exists well-posed s-periodic initial value problems, for almost all 

Ph ' s £ iP. From the Lax-pairs we furthermore derive invariants for corresponding 

reductions to dynamical mappings for some explicit examples. 
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^ ■ 1. Introduction 

'nI ' The discrete-time Toda equation, i.e. the time-discretized version of the usual 

^~. , Toda chain, which is given by the foUowing partial difference equation (PAE) 

in _ 

^P , Tl~l,m~lTl+l^m + l = Tl + l.rn-lTl-l.m+l ~ Tl ,m (1-1) 

O 

plays an important role in many areas of mathematical physics. It is probably the 
first integrable fully discrete equation that can be found in the literature: it ap- 
i^j ■ pears, albeit in slightly different form, for the first time in the paper by Frobenius 

i*^ , [13] as an identity for certain determinants used in the determination of rational 

C^ ' approximations of functions given by power series. The latter are nowadays known 

as Fade approximants, which, in fact, were introduced by Frobenius 10 years prior 
to Fade's work, (cf. [M] for a more modern account). The discrete-time Toda 
equation re-entered the modern era through the literature of integrable systems 
|18j . Hirota in |18j introduced the equation as a natural time-discretization of the 
famous Toda chain equation, generalizing the latter to a partial difference equation 
on the space-time lattice (here the lattice sites are labelled by the discrete indepen- 
dent variables (^tti) G Z^). Eq. (jl.l[) exhibits the prominent integrability features, 
such as the existence of multi-soliton type solutions and the existence of an under- 
lying linear problem (Lax pair). This FAE has also appeared in physics, namely 
as the nonlinear equation governing the spin-spin correlation functions of the two- 
dimensional Ising model, cf. [26[ I31j . The connection we investigate and generalize 
in this paper, is the emergence of (jl.ip in the theory of formal orthogonal polyno- 
mials, where the equation ()1.1|) is connected to the well known quotient-difference 
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(QD) algorithm of Rutishauser, [33], cf.[30]. 

ei,m+l + ft+l.m+l = <ll+l,m + ei+l,Tm 

ei,m+iqi,m+l = qi+l,mei^m- (1-2) 

The theory of formal orthogonal polynomials is a subject central to modern 
numerical analysis, in which orthogonalities are investigated on their general prop- 
erties regarding the recurrence structure from a formal point of view (i.e., without 
specifying particular classes of weight functions), with a sight on the construction 
of numerical algorithms, rather that on the analytic properties (such as the moment 
problem, or the problem of the behaviour of the zeroes of the polynomial) arising 
from the particular properties of the weight functions. In this area of research gen- 
eral constructions such as those of vector Pade approximants, adjacent orthogonal 
families (where connected sequences of orthogonality functionals are postulated), 
and associated convergence acceleration algorithms, and factorisation methods have 
been developed, cf. e.g. [SJ[T2]. In this context, the QD algorithm emerges as a 
prominent method to locate the zeroes of analytic functions or to compute con- 
vergence factors for asymptotic expansions through formal power series, using a 
finite-difference scheme and continued fractions, cf. |17j . It was pointed out in |30| 
that both the QD algorithm and certain convergence algorithms, |15| . are intimately 
related to integrable discrete systems. In fact, the famous £-algorithm of Wynn, 
[41] , surprisingly turns out to be identical to the lattice potential Korteweg-de Vries 
(KdV) equation, a well-known exactly integrable PAE, thus allowing us to inter- 
pret this numerical algorithm as a symplectic dynamical system with a rich solution 
structure. Siinilarly, the famous "missing identity of Frobenius" appearing as the 
rhombus rule in the Pade tables and governing the stability of the e-algorithm, cf. 
[32], can be identified as an exactly solvable PAE closely related to the Toda lattice 
and discretisations of the KdV equation. 

In this paper we introduce a novel integrable PAE: 

f/-(-l,m-2f/_l,m+l {(Jl,m+2'^l,7n-l — <^ I, mO' I, ni+l) 
= O'l^m-lO'l-l^m+2 i<^l-l.m+lO'l+2.m-2 — 0'Lm'^l + l,rn~l) (1-3) 

+ 0'l^m+l<^l + l,m-l {o'l,m-2'^l-l,m+2 — <^l + l,7n-2<^l-2,m+2) , 

which we argue can be regarded as a Higher order Analogue of the Discrete-time 
Toda equation and therefore name it the HADT equation. It is defined on a stencil 
of 11 points in the lattice as depicted in Figure [TJ 



Figure 1 . The stencil of the HADT equation ([L3]) 

The derivation of the HADT equation is parallel to the way in which the discrete- 
time Toda equation emerges in the theory of orthogonal polynomials: in fact, we 
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introduce a family of two-variable orthogonal polynomials restricted by an elliptic 
curve and exploit the recurrence structure for the relevant Hankel determinants. 
The general problem of orthogonal polynomials on algebraic curves was discussed 
in the monograph by Suetin, |351 Chapter 7], and has also been addressed in the 
context of the study of formal orthogonalities, cf. e.g. j8]. The construction 
of the associated two- variable orthogonal polynomials on elliptic curve is pursued 
in |34| . whereas in this paper we will concentrate on the HADT equation ()1.3|) 
itself and its reductions. Nevertheless, we will present the derivation from the 
adjacency structure of elliptic two-variable orthogonal polynomials, as it produces 
not only the nonlinear equation itself but also its Lax pair. Furthermore, the 
recurrence structure of the relevant Hankel determinants also yields a novel type of 
QD formalism, namely the following so-called quotient-quotient-difference (QQD) 
scheme: 

Mi+2,m + I'i+l.m + W(+l,m+l = Ui^„i+3 + U;+i^m+i + Wj+i^m, 

Ul,m+3Vl.m+l = Vl+i^r,iUl+l,m, (1.4) 

which we believe may have future applications in numerical analysis. 

Once we have obtained the HADT equation, the QQD scheme and their Lax 
pairs, we address the issue of periodic reductions. A periodic solution of a (system 
of) lattice equation(s) is a solution which satisfies Ui.,„ = Ui-^si.m+s2 for some 
si,S2 G Z. By imposing periodicity a lattice equation reduces to a system of 
ordinary difference equations, or, equivalently, a mapping. The so-called staircase 
method utilizes the Lax-pair of the lattice equation to construct integrals for the 
mapping. The method was introduced in [29], where it was applied to the KdV 
equation and to a mixed modified KdV-Toda equation, taking si = S2 and si = 
3-2 + 1. In [32] a two-parameter family of periodic reductions was studied, with si 
and 32 being co-prime integers. Recently, in [37] a unified and geometric picture 
for periodic reductions, with nonzero s = (si,S2) G Z^ has been provided. In 
[37] it was shown, for any given scalar lattice equation on some arbitrary stencil 
of lattice points, there exists a well-posed, or nearly well-posed, s-periodic initial 
value problem, for all nonzero s G Z^. Therefore, the trace of the monodromy 
matrix provides integrals for any periodic reduction of any intcgrable scalar lattice 
equation, and it is expected that the same is true for systems of lattice equations, see 
[37] [38] . In this paper we generalize the approach of [37] and develop a systematic 
method to construct well-posed periodic initial value problems for systems of lattice 
equations. This method is applied to the QQD scheme as well as to an intermediate 
system of PAEs, cf. equation p.l9p . 

On general principles, one expects periodic reductions of integrable lattice equa- 
tions to lead to integrable maps in a precise sense, namely to complete integrability 
in the sense of Liouville- Arnold [9] [39]. This requires, in addition to the existence 
of a sufficient number of integrals, some global properties as well the existence of 
a symplectic structure which is preserved by the discrete dynamics, and with re- 
spect to which the integrals are in involution. The latter are issues which we are 
not addressing in the current paper, where we have developed the basic structures 
and classified the consistent reductions from the lattice equation together with the 
explicit form of the integrals. 



4 PAUL E. SPICER, FRANK W. NIJHOFF, AND PETER H. VAN DER KAMP 

Furthermore, there is the problem of counting the number of independent inte- 
grals to ensure that in principle a 1-1 map to the relevant number of degrees of free- 
doms exist for the mappings concerned. In fact, typically some reduced variables, 
i.e. suitable combinations of the variables on the vertices, need to be introduced, 
leading to a dimensional reduction of the mappings. This issue is addressed in the 
section 5. For the reductions we obtained from the HADT equation, as well as those 
from the QQD scheme, we are able to find a sufficient number of independent inte- 
grals, that is, equal to at least half the dimensionality (i.e. number of components) 
of the mapping. For the reductions obtained from the intermediate system p.l9p 
this, however, is not the case. In j38| a systematic method was presented, which 
exploits symmetries of the (system of) lattice equation(s) to reduce the dimension 
of the mappings. For all cases that were considered in [38], the staircase method 
provides a sufficient number of integrals for the dimensionally reduced mappings to 
be in principle completely integrable. We arrive at a similar conclusion, however, 
in certain cases the staircase method provides a product of 2-integrals J J' ^ from 
which a missing integral can be obtainedu 



2. The discrete-time Toda equation and the QD algorithm from 

orthogonal polynomials 

We present here the derivation of the bilinear discrete-time Toda equation us- 
ing formal orthogonal polynomials. Recurrence relations for adjacent one- variable 
orthogonal polynomials provide a Lax-pair for this equation, as well as for the as- 
sociated QD algorithm. The results presented here are known, see for example j30j , 
however they outline the main idea which we employ in the setting of two- variable 
elliptic orthogonal polynomials in the next section. 

Let P„ = Pn{x) be a family of polynomials of a variable x orthogonal with respect 
to a weight function w{x) on a curve F in the complex plane. The weight function 
defines a linear functional L with respect to which we can define the moments c„ 
as 



£(e„ 



x^w{x)dx i = 0, 1, . . 



(2.1) 



where e„ are the monomials. The P„'s, which are assumed to be monic, satisfy a 
three term recurrence relation of the form 



^^-i n — t^n+1 I i>nJ^^n i -flji-fVi— 1 



(2.2) 



(where Sn and i?„ are finite) and in order to insure their existence we assume that 
the corresponding Hankel determinants are non-zero. 

We introduce a family of adjacent orthogonal polynomials (which are orthogonal 
with respect to a shifted weight x™'w{x)) [3]: 



p(r\a 



Am) 



^n+m — l 

1 



Cn+r, 



C2n+n 



(2.3a) 



Recall, a function J is an fc-integral, or fc-symmetry, of a mapping if it is an integral, or 
symmetry, of the fc*'' power of that mapping |16| . If one has one fc-integral, then one can construct 
k of them, or, even better, k integrals. For example, it is easy to see that J" = J implies that 
both J J' and J + J' are integrals. 
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with the corresponding Hankel determinant: 



and the two row/column Sylvester identit}o [4l[27], 



^-■n-j-m 



(2.3b) 



(2.4) 



Applying two different forms of the Sylvester identity to the determinant for P, 
leads to the derivation of two xPn recurrence relations 



(m) 



P 



p: 



(™) ^ ^p(™+l) „ T/(™)p(™) 
a+l n *^ n -^ n ^ 

(™) - ^p(™+2) _Tl/(m)p(m+l) 



)l + l 



with 



]/("*) 



(m+l) A (m) 
. (m+1) . (rn) ' 



w(™) 



A, 



(m+l) A (m+l) 



A (m) A (m+2) 



(2.5a) 
(2.5b) 

(2.6) 



The combination of these relations leads to the nionic recurrence relation, of the 
form (12:21). 



From this approach the coefficients S'„ and i?„ can be further simplified (in terms 
of the Hankel determinants): 



On — 



-Rn = 



A^ 



(m+l) . (in) 



(m) A (m+l) 



f.(m).(m) 
^n ^n-2 

A (m) A (m) ' 
^n-l^n-1 



A (m+l) . (m) 

A (m+l) A (m) ■ 
^n-1 ^n-1 



(2.7a) 



(2.7b) 



where we have suppressed the m-dependence in the symbols i?„ and 5„. We achieve 
this simplification by making use of a bilinear relation that exists between the 
Hankel determinants A„ . This bilinear relation is found by applying the Sylvester 
identity to A„ ^ 




A(m) A(m+2) _ A (m+2) A (m) a (m+l) a (m+l) 

^n ^n-2 ~ ^n-1 ^n-1 ^n-1 ^n-1 



(2.8) 



This identity has many different names including the Jacobi identity, Lewis Carroll's identity 
and the window-pane identity, however we will just refer to it as the Sylvester identity. 
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We will refer to equation (|2.8p as the discrete-time Toda equation. Indeed, 
substituting A„ — r_2,i-m-i,m-i S'lid taking I = —2n — to as new independent 
variable yields the more standard form of this equation (|1.1[) . The continuum limit 
leading to the usual Toda chain is most easily seen starting from equation (jl.ip . 
multiplying with ipq and applying a point transformation ti ^ '-^ Ti^mCt^ /3'™ with 




/3 = 

leads to 

(P - q)'^Ti-i^m-iTi+i,m+i - {p + 9)^T;_i,,„+ir;+i,„j_i + ^pqrf,^ = 0, (2.9) 

in which p, q are lattice parameters (associated with shifts on the lattice in I- and in 
directions respectively). Performing a changes of variables, {l,m) ^^ {n ^ 1 + m, m) 
and taking a limit I — >■ oo, m — >■ oo, n fixed, g — p = e — > 0, em -^ t, Taylor 
expansion yields n^m = T„(io + e"i) => t;,„,+i = t„+i + e-r„+i + . . . and we get 
the semi-discrete bilinear equation: 

which is related to the usual Toda chain equation by setting (f>„ = log(T„/T„_2)- 

Thus the shadows of integrability already appear in the underlying structure of 
the standard theory of orthogonal polynomials. 

2.1. A Lelx pair for the QD algorithm and the discrete-time Toda equa- 
tion. We now view equation (|2.8p as an intcgrable lattice equation, as opposed to 
an identity for Hankel determinants (j2.3b[) . To this end we write A„ -^ = t„ 



' n.7n-j 



and similarly Vn"^' = Vn„i, Wn = Wn^m- The variable x will play the role of 
spectral parameter and will be denoted x ~ X. 

The relations (|2.5p constitute a Lax-pair for equation (j2.8|) . Working with 
the fields v,w we first derive a related quotient-difference system. Let ^„_m = 
{Pn , Pn ); then using the recurrence relations we derive 



ith 



^n+l,m — Ln^m^n,m ^nd ^n,m + l — -A'^n,m^n,m, 






(2.10a) 



Mn,m = ( .X 1 + A-l^ )■ (2.10b) 

The compatibility of the two linear systems is then equivalent to the discrete Lax 
equation: 

= Ln^m+lMn^m ~ Mn+l.mLn.m 



-X^^Vn,mE E + X^'^{Wn,m{Vn+l,m + E) - Vn,m+lWn+l,m) 

where 

E = Vn,m+2 — Vn+l,m + W7i+l,m — Wn^m+1- 
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Thus L,M provide a Lax pair for the quotient-difference scheme 

Vn,ra+2+ Wn+l,m = Vn+l.m + Wn^m+li (2.11a) 

Wn,mVn+l,'m = Vn.rn+lWn+l,m- (2.11b) 

This scheme is related to the (more standard) QD algorithm (jl.2p . by qn,m ~ Vn.m, 

Upon substitution of p.6p into (|2.11ap this equation is a consequence of p.Sp . 
whereas substitution of (|2.6p into (|2.11bp turns into an identity. To find a good Lax- 
pair for equation p.Sp . one substitutes equations p.6p into the above Lax matrices, 
and one uses (|2.8p to simplify 

^n—l,rn-\-l'^n~\-l,m 

Using the same simplification, equation p.llap can be written 

= En.m+l + Wn+l,m ~ Vn+l,m 

_ Tn+l,m+l / _ / s2\ 

'^n,m-t-2'^n,T7i+l'^n + l,m 

where one recognizes the discrete-time Toda equation ()2.8p . 
3. The HADT equation and a QQD scheme from elliptic orthogonal 

POLYNOMIALS 

In this section we derive the HADT equation. We follow a similar route as in the 
previous section, except that here we consider two- variable orthogonal polynomials, 
where the variables are restricted by the condition that they form the coordinates of 
an elliptic curve. The recurrence relations for these elliptic orthogonal polynomials 
yield both a Lax pair and the QQD scheme. 

3.1. Two- variable elliptic orthogonal polynomials. As a starting point for our 
construction we introduce the sequence of elementary monomials |34| . associated 
with a class of two variable orthogonal polynomials. We consider a sequence of 
monomials where the x and y are given different weight, namely x ^ 2 and y ^ 3: 

eo = 1, 62 = a;, ea = y, 64 = x^, e^ = xy, eg^x^, ■■■ 

or, in general: 

eo{x,y)^l, e2k{x,y) = x'', e2k+i{x,y) ^ x''~''^y, A: = 1,2,.... 

In comparison with the polynomials of Krall and SchefFer |21| , their two- variable 
orthogonal polynomials consist of x '^ 1 and y ^ 1 oi equal weight. The reason for 
choosing different weights is the variables x and y are related through a Weierstrass 
elliptic curve. We use this sequence as our basis for the expansion of a new class of 
two- variable orthogonal polynomials taking the form: 



Pfc(x,y)=^pfe,(x,2/), (3.1) 



(k) 

which arc monic if the leading coefficient pj^ = 1. 
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In the spirit of the formal approach to orthogonal polynomials, of. e.g. [3 [5], we 
assume that a bilinear form (, ) exists and can be derived from a linear functional 
L and consequently we can define the associated moments by 



Cfe = £(efc) 



(3.2) 



The assumption of the existence of an inner product (, ) on the space V spanned by 
the monomials e^, is such that 

{xP, Q) = (P, xQ) 

for any two elements P,Q & V. 

Under the assumption of orthogonality and using the standard Gram-Schmidt 
orthogonalisation (through the use of Cramer's rule), leads to the following expres- 
sion for the adjacent elliptic orthogonal polynomials: 



H'\^,y) 



(e;,eo) (e;,e2) 

(ei+i,eo) (e;+i,e2) 



(e;+fc-2,eo) (e(+A;_2,e2) 
eo 62 



(e(,efc) 
(e(+i,efe) 



efc 



together with the corresponding Hankel determinant: 



M'ii , 1^0 



(3.3) 





{ei,ea) 


{ei,e2) 




{ei+i,eo) 


(ei+i,e2) 


Af = 








(e;+fc_i,eo) 


(ei+fc_i,e2) 


gre for / 


= 0: 





P!.'\^,y) 



(eo,eo) (60,62) 

(62,60) (62,62) 



(6A:-1,6o) (6A;_i,62) 
60 62 



with 



A 



(0) 



(60,60) (60,62) 
(62,60) (62,62) 



(6fc,6o) (6A;,62) 



{ei,ek) 

(6i+i,6fc) 



(6i+fc-l,6fc) 



(60,6fe) 
(62,6fe) 



(6fc-l,6A;) 
6fc 



(60,6fc) 
(62,efc) 



{ek,ek) 



1^0 , 



(3.4) 



/^t'\ , (3.5) 



(3.6) 
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. In addition to the polynomials p.3p . we also need to introduce the polynomials: 



i(0 



'k'i^^y) 



(ei,eo) (e/,e2) 

(ei+2,eo) (e;+2,e2) 



(e;+fc-i,eo) (e;+fe_i,e2) 
eo 62 



together with its corresponding Hankel determinant: 



{ei+2,ek) 



{ei+k-i,ek} 
e-k 



KU 



(3.7) 



e 



(i) 



(e,,eo) (e/,e2) 

(e;+2,eo) (e/+2,e2) 



(e;+fc,eo) (e/+fe,e2) 



(e/+2,efe) 



(ez+fc,efe) 



(3.8) 



noting that 



gf = n^"' 



ei°) - Af 



Remark: We note that for / ^ 0, 1 the polynomials Pf, do noi form an orthogonal 
family. In fact, from the determinantal definition p.3p we immediately observe that 



(e,,p(") = (e,+i,p(')) 



(e,+fc_2,P^")=0 , {ei+k-uPf^) 



.(0 



^(0 



whereas 






A 



(0 
fe-i 



Wc now proceed using determinantal identities of Sylvester type (Appendix A) 
to derive relations between the polynomials P^: and the Hankel determinants. 



3.2. Recurrence relations in PJ: and Qj^ ■ Using a 3 row/column Sylvester 
identity wc find it is possible to find x-rccurrence relations and linear relations in 
Pj. and Q\. ■ To achieve the former it is necessary to fix the columns, so that 
eo = 1 and e^ ~ y are removed from the determinant and the position of the 
row removal is dependent on restricting the introduction of new objects. Hence we 
apply the following cutting of three rows and columns (|A.6p to the determinant for 

n(0. 
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(where the small red line indicates a space between the first column and the third) . 
We are led to the recurrence relation 



^k ~ ■^^k-2 ^k-2^k-l 



Wt2Pk 



C + l) 
1 



1^0,1, 



(3.9) 



where 



For 



V, 



(0 



we have 



P, 



(0) 



A (I) A (1+3) 

^k ^fc 

a(0 Ail+3)-' 

^fc+l^fc-1 



^, ^ Ai!-)Ar^ 



xP, 



(4) 
fe-2 



^(0) ^(4) 



P, 



(0) 



fc-1 



A 



(0 A (1+3) 

-l^fe-1 



^(2) 0(2) 
^fc-2'^fc-2 
^(0) ^(4) 
^fc-l^fc-3 



P 



(2) 



fe-1' 



(3.10) 



(3.11) 



whilst obviously, since P^ is not defined, there is no relation for Z = 1. 

By making use of intermediate determinant expressions (Appendix [B]) we can 



derive additional relations in terms of P, 



(0 



and gf . 



Thus in addition to the 

(0 



recurrence relation p.9p we have also derived the following relation in terms of P^ 
feauation lB.ip : 



3(0 



P, 



e+i) 



[/, 



(0 



in 



(;+i) 
1 ' 



with 



u: 



(0 



aW aC+i) 



Similarly, we derived recurrence formulae for the set {P^ 
equations (|B.9p . ()B.lip and (JB.ISP respectively), 



W 



P. 



xQj^_2 



aW C)('+2) 
^fc-2^fc-2 p(0 
(l) o(i+2) -^fc-] 



e 



'2A, 



A^iie 



A;-3 



aC) «('+2) 
^A;-l*^fe-3 



?i'^ 



_ p(i+l 



P'' 



1) -^fe-l 



QL' 



w 



p 



(0 



a(0 a(' + 2 

0(0 ^{i- 



A(O0(O 



('+2) 



^fc-l^fc-l 



(3.12) 



(3.13) 



Q^l^} (Appendix 1131 



(1+2) 



"^fe-l' 



/^0,1 



(3.14) 
(3.15) 

(3.16) 
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3.3. Bilinear Hankel identities. In appendix IB. 21 we construct a pair of bilin- 
ear relations: a four term bilinear relation (derived by applying a 3 row/column 
Sylvester identity to the determinant Aj, ), and a three term bilinear relation (re- 
sulting from the combination of other Hankel identities) : 



^^^^ = 


= A«,Ai?^) 


^fe-l^fc-2 +^fc-l ^k~2 ' 

(3.17a) 


Ai'^Ai?-) 


= ei?'Ai?-)- 


-e'^l^A'i^^l (3.17b) 



3.4. The HADT equation. The system p.l7p can be rearranged to eliminate 



the 6j, . Thus we have 



a(0a('+4) .(0 .(;+4) „(i+2) 
^k ^fc-3 ^fe-l^fc-2 *-'fc-2 


*-'fe-i 


.{l+l).{l+3) A(i+l)A(i+3) .{1+3) 
^fc-1 ^fe-2 ^fc-1 ^fe-2 '^k-2 


^fc-i 


Ai?'Ai^f) ei?) 


*-'fe-i 


.(/+l)^(i+l) A(i + 1) 


A C + l) ' 
^fc-l 


which can be expressed in a simpler way using 




xH' = T^-rt. 




^{l) _ -p(i+2) p(0 
^k — '' k-2 ^fc-1' 




where T^'^ = e^'VAfc"*"^^. From these two expressions we 


have 




p(0 


yil+2) I vC) ^ yC) - r('+2) 
^fc+l ' -^fe+2 ' ^fe+1 ^ -^ fe+1 




and by rearranging, we find a quadrilinear relation in terms of A^^ on an 11-point 


stencil, the HADT equation: 




Ai'jj (A('+^)A('+^)Ai'+,^) - Ai'+2)A('+^)A[;+f 


I a('+4) aO+S) .(i-Hl)\ 
+ ^fe-2 ^fc ^fe+1 J 


= Ai?^f(Ai^,Ar^'Ai?-)-Ar)Ar)Al?:^) 


+ Ai?)Ar)Ai^)) . 



(3.18) 

The HADT equation bears its name as a Higher order Analogue of the Discrete- 
time Toda equation, due to the similarities between how both equations arise in 
the theory of formal orthogonal polynomials. We expect the HADT equation to 
also be integrable and derive a Lax pair for it in the next section. Also note the 
similarity between the stencils of the two equations (j2.8p and p.lSp . as displayed 
in Figure [21 

It would be interesting to know how to enter lattice parameters (i.e., parameters 
associated with the grid steps) into the HADT equation, as we did for the discrete- 
time Toda equation, see equation p.9p . and then to study its continuum limit. 
However, there is a lot of freedom here to enter lattice parameters through point 
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•-.■ • -f 

•■ # • 



* -• 



Figure 2. The stencils of the discrete-tnTie Toda equation 
and the HADT equation ((3?T8| 



transformations, and this should be properly investigated. We leave that issue to 
a future work. 

Next we follow a similar approach as in section 12.11 and derive Lax pairs for 
the several Hankel identities viewed as integrable systems, using the recurrence 
relations for the two- variable elliptic orthogonal polynomials. As before the fields 



will carry two sub-indices; we denote x = A, Aj™ = (Ji.m: 6; 



Pl,m, U^. 



(0 



Uk,h 



etc. By not evaluating the coefficients in the recurrences ()3.12p and (|3.9p . we will 
find a novel quotient-quotient-difference scheme. 



3.5. A Lax pair for the system of lattice equations ()3.17p . We derive a Lax 
pair for the system of equations p.l7p which in terms of cr, p reads 



0'l+l,m-2<7l-2,m+2 + Pl.m,-2<7l-l,m+l 



0'l,m-lO'l-ljn+l + Pl-l,m-lO'l,m 



0'(,m-20"i-l.m+2 + <^Lrn-lPl-l,m, 



Pl,'^ 



-lO'l-l.r, 



(3.19a) 
(3.19b) 



Taking *,,„ = (Q|'"+'\p/;;+'\g[;"i\p/"2^) we derive *z+i.™ = i^™*,,™, and 

^l+l,m = Mi^rn^i^rn, with 



Li, 



+ X 



^l + l,m + lPl-l,m + l 
CTi,m+lPi,m+l 










^l + 1.7Ti + 2Pl-l,m+l 



— ^l + 2,m^l.r 


n + 1 


0"l + l,mCTi + l, 


ni + 1 







— ^l + l,Tn + 2'^l 


m + 1 







— °'i + 2,mPl,j 









Pl + l,m + lPl,r, 



^I,m + 3'^/ + 2,m 







-Pl + l,m + lPl,n 
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and 



/ —pi,m + 2(^l,m + 2Pl-l,m+l 



M, 



l.m 








i + lPi-l,m+l 





-Pl.m + 2Pl,m + lPl,m 

-l,m + 2f^l + l,m(^L,m + 3 
Pl + l,m<yi,m. + l 

<^l + l,m<^l,Tn + 2 
1 



Pl,m + 2Pl,m+lPl,r, 



Pl- 



— gj + l.iTi+lPi.m 
<^l + l,m.<7l,m.+2 





+ A- 






Vo 



o-|, 












1/ 



1 + 2 CT/ . + 1 . 











Pf ,m + 2C"Z-l,m. + 3 







/ 



Here we have used the equations p.l9p to reduce the number of terms in the entries 
of the Lax matrices as much as possible. The compatibiUty condition Li^m+iMi^m — 
Mi^i,,nLi.,n = is equivalent to the system (|3.19p . However, to verify that the left 
hand side vanishes modulo the system, some work has to be done. The equations 
do not factor out as nicely as in section [Ol 



3.6. A Lax pair for the HADT equation, and the QQD scheme. We first 

derive a Lax-pair for the following QQD scheme, satisfied by the coefficients (C/^ " ~ 
ui^rn etc.) in the recurrence equations (|3.9p and (|3.12p : 

Ul+2,m + Vi + l,rn + Wl+i^rn+l = Ui^rn+Z + Vl+l,m+l + Wl+l,rn, (3.20a) 

Wi,m-f3Wi,m+l = Vl + i^rnUl + l^m, (3.20b) 

ULm+2Wl^m = Wl+l,mUl+l^m- (3.20c) 

Let*; '"""' '""'"' '""'" '~~' 



'u-a = {Pr \Pi+i \Pi+2 ' ^Pi+i)- We use ^M and ^M to find 
*i+i,m = Li^m^i^m, and *i,,„-fi = Mi^rn'^i^m, where L and M are given 



V 



-Ul.m+2 







1 






/ 



+ A 



and 



\ U/-|-l,m-|-2UZ,m-|-2 



Ml, 



\ 



Ml,Tn + 2 

-Ul^m+2 






/ -W;,„: 





V 





1 

— Ul+2,m 
— Wi+2,m1J-l+2,r 





M; + l,m-|-l — Ui^i, 





1 






1 

-Vl+2,rn + Wl+2,: 

\ 















1 


J 



1 



-ui+i, 
























Here we have used the scheme (|3.20p to simplify the A~^-part of Mi,m- The com- 
patibility of these linear systems is equivalent to the QQD scheme p.20p . Note, by 
substituting the expressions p.l3p . p.lOp for u,v,w in terms of a into the QQD 
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scheme, the first equation is equivalent to the HADT equation, whereas the later 
two are satisfied identically. Therefore, after the same substitution, the matrices 
L, M provide a good Lax pair for the HADT equation p.lSp . The stencils of the 
Lax matrices L(A),M(A) are depicted in Figures |3l S] 



. ^.™ . '■■. ::* 

Figure 3. 

The stencil of L;,™ Figure 4. 

The stencil of M/^m 

4. Periodic reduction 

In |37| a method was given to obtain initial value problems for scalar equations 
on arbitrary stencils. The question was raised whether a similar construction can 
be done for systems, and one example was given (the quotient-difference algorithm) 
where this is the case indeed. Here we will present a systematic approach towards 
constructing initial value problems for systems of partial difference equations. 

Using the method of [37] it is easily shown that the dimension of the s-periodic 
initial value problem for the HADT equation is given by the piecewise linear func- 
tion 

4max{|si+S2|,|si|}. (4.1) 

We will show that the s-pcriodic straight band initial value problems for the in- 
termediate system p.l9|) as well as those for the QQD scheme are of the same 
dimension as the initial value problems for the HADT equation, namely (|4.1|) . 

Although we do give a method to construct s-periodic initial value problems for 
systems, the question remains whether this can be done algorithmically. As we will 
see, one can not easily determine the dimensionality of the reductions, nor which 
regions for s should be distinguished, by simply looking at the stencils on which 
the equations of the system are defined. 

4.1. Initial value problems for the HADT equation. To pose initial value 
problems for equation (jl.3p we use the method developed in |37| . By the S -polygon 
of the equation we mean the boundary of the convex hull of the stencil S of the 
equation (the dotted parallelogram in Figure[5]). The S -directions are the directions 
in the S'-polygon of the equation: (1,0) and (1, —2). Therefore, we distinguish two 
regions in the plane 

R^ = {seZ^ : a> 0,6> 0,e = 1 or 6> 2a,e = -1}, 
i?2 = {s eZ^ : < 6< 2a,e= -1}, 

as in Figure [S] If s S i?^ or s G i?^, then performing an s-periodic reduction yields a 
well-posed initial value problem, i.e. the s-periodic solution is obtained by iteration 
of a single valued mapping. If on the other hand s is on a boundary between R^ 
and R^ one needs to impose periodicity on the solution, and then one gets a multi- 
valued mapping, or correspondence. To each of the regions i?' we associate a vector 
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Figure 5. The 11-point stencil of the HADT equation, its convex 
hull, and two regions in the plane 



d' which is a difference between the two points of 5*, where a line with direction 
s e i?' intersects the S'-polygon in a single point. For the HADT equation we have 

di = (0,4), d2^(4,-4). (4.2) 

Hence the dimension of the s-periodic initial value problem is, cf. [37l Equation 8], 



maxj=i,2{|Dct( 1|} = 4max{|si + S2|,|si|}. 

We give two examples before we present the general s-periodic initial value problem 
in terms of s-reduced variables. 

4.1.1. (2,- l)-reduction of the HADT equation. We specify initial values between 
two parallel lines that squeeze the stencil as in Figure [HI Clearly, given the 8 

,a;3. .,X2. .^1 




s^^g ,xs 



X7 -.xe X5 



Figure 6. (2,-l)-periodic initial values 

values xi,X2, ■ ■ ■ ,xs at the black dots one is able to determine the value xg at the 
white dots on the fat-dotted line, using the HADT equation. Thus we get an eight 
dimensional mapping: 

Xi H- Xi+i, 1 < i < 8, 



Xs H- 



a;ia;4a;7X8 x^xj x^xg X4a;5a;8 x^x^ 



(4.3) 



X2a;3a;6 ^3 X2 X2X& x^Xf, 

We can calculate integrals for this mapping by constructing a so called mon- 
odromy matrix, which is a product of Lax matrices along a staircase over a one 
period long distance. Taking the point at the bottom left corner to be the origin 
{I, m) = (0,0), we define the monodromy matrix to be the inversely ordered product 
of Lax matrices along a staircase from (/, m) = (1, 0) to {I, m) = (3, —1), 

L = A/3 _j^i2, 0-^1,0- 
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Using the stencils of the Lax matrices, see Figures [3] and IH the reader can verify 
that the matrices in £ depend on initial values only. The characteristic polynomial 
of L is 



Det(AJ -L)= ^i^i^i - Xf + 1) + /i/i(/i2 -Xfi-X)- A(/2M^ + /3/i + 1), 



where 
h 



xsxe , X2X^ xixs, 



X1X4 



XsX5 X^Xi 



X2X7 



Xt 



XqXz 






XqX^ 
XiX-j 



X2X3 XqXi X2X7 



a;3X5 



XiXQ 

X2X5 , 



X2X4 

XqXs ^ 



X2XG 

X^Xy 



)(- 



xt 



Xa 



X2Xii 



X2XgJ \X5X7 

xixs X2X7 



X3X5 a;3a;7 



X4Xe\ 
X3X7J ' 



a;3a;4 a;4a;5 x^xg xaXe xsXq 

Because we subtracted /i from the coefficient of /i^A in the definition of I2 it 
factorizes nicely as J J', where J is a 2-intcgral of the map and J' is the image of J 
under the map. Remember a fc-integral is an integral of the fcth power of the map 
[16]. Indeed, given that J" = J, we find that the product J J' = J" J' = {J J')' 
is an invariant. But now, we can write down another (functionally independent) 
invariant, namely J + J' . Moreover, the integrals Ii, I2, I3, li, with 



h 



+ 



X1X7 2:3X5 



X22:8 XiXe 



X4Xe X2Xi X2Xe 
are functionally independent. 



X2Xe x^xj X3X5 xsxj 



xzx^ 



4.1.2. (2,-3)-reduction of the HADT equation. We specify initial values between 
two parallel lines that squeeze the stencil as in Figure [T] 



J?3..J>?i5..j:9 




,X3 -^xe -.^xg 
^xi_^X4 -^xr 
Figure 7. (2,-3)-pcriodic initial values 



Clearly, given the 8 values xi , a;2 , . . . , a;8 at the black dots one is able to determine 
the value Xg at the white dots on the fat-dotted line, using the HADT equation. 
Thus we get an eight dimensional mapping: 

1< i < 8, 



X-j ^~r X 



i+li 



xs 1-^ 



Xixr 



X5 



X7X5 
X3 



(4.4) 



1-4 V 2:3 

To be able to evaluate monodromy matrix 

we first have to determine the value xio- Using the stencils of the Lax matrices, 
see Figures [3] and SI it seems we also need to to determine the value xo , as both 
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io,i and Mf^i depend on xq. However, their product does not depend on xq. The 
characteristic polynomial is given by 



Dct{^ll -L)= fi^ii^ - A)2 + (m - A)/2 - ^iJJ' + 1) + (/i - A)''' - /i(M - A)2/. 



where 

h 



J 



X4X5 X3XS XiXq X2X7 

XiX& X^Xq X3X4 XiX^ ' 

3:23:7 ( xxx\ x\xfi , -,\ , 3;iX8 _ X6a:2 _ x^x-j 



XzX(, \ 2:4X5 X4X\ ) x^x^ x\ x\ 

a;ia;7 - X5X3 



where J is a 2-integral of the mapping. The four integrals Ii,l2, J J' and J + J' 
are functionally independent. 

4.1.3. General s-reduction. By imposing s-periodicity, that is, 

our PAE reduces to a system of r — gcd(si,S2) OAEs. We perform so-called 
s-rcduction. which is basically a change of variables {I, m) — > {n,p). 

Definition 1. We define a,b,c,d G N and e = ±1 by a/b = |si/s2| = esi/s2; 
gcd{a, 6) = 1; 6 = ^ c = 0, d = 1; 6 = 1 => c = 1, d = 0; fe > 1 => c is the smallest 
positive integer such that be — ad ~ 1 . Now we specify 

ni^m := bl — earn, Pi,m '-^ ^dl + ecm mod r, (4-5) 

The variable n tells us on which line with direction (a, eb) the point is and p 
distinguishes the r inequivalent points on each line, see also [37]. The important 
properties are 

("-,P)i+a,rn+£6 = {n,p + 1);^,„, ("-, p)i+c,rn+£d = ("- + l,p)z,,„ 

where the latter shift is the standard way to update the initial values. 

We adopt the notation of [37] , where reduced variables carry one upper and one 
lower index, that is, the reduction will be denoted ai^m ^^ <^n- When specifying 
(periodic) initial value problems in terms of reduced variables we specify sets of 
variables a^ where n runs over a specific range and p runs over the full range 
pGNr :={0,l,...,r-l}, cf. [37]. 

The range of n is determined by the differences associated to each region. For the 
HADT equation, when s <E R^ (where e = —1, see Figure |S]) the range of n will be 
fromO to n{d^)-l = fe-O + a-4-1, and when s £ R"^ we take < 71 < 4(6 + ea)-l, 
using the differences (|4.2p . 

Thus we obtain 

Proposition 2. The HADT equation admits a well-posed s-periodic initial value 
problem if (a, 6) is not equal to (1,0) or (1, —2). 

• For s ^ R^ the set {a^ : n £ ^4aiP G ^r} provides a well-posed initial value 
problem of dimension 4| si |. 

• For s G i?^ the set {aP : n G N4(h_|_(:a), p G N,} provides a well-posed initial 
value problem of dimension 4| si + S2 |. 



18 



PAUL E. SPICER, FRANK W. NIJHOFF, AND PETER H. VAN DER KAMP 



Remark: In section H.l. li the initial values {ctq, . . . , cr"} are conveniently denoted 
(tJJ = xs-n and the mapping is the backward shift n i— > ?i — 1, or (Z, m) t-^ (^ — 1, rn). 
In section r4.1.2l the initial values {ctq, . . . , cry} arc denoted cr° = Xn+i and the 
mapping is the standard one n ^^ n + I, that is (Z, m) t-^ (I + l,m — 1). 



4.2. Initial value problems for the system of lattice equations ()3.19p . The 

stencils of the system p.l9p are depicted in Figure El The S'-directions in the 
^-polygons of the equations are (0,1), (1,0), (1,-1), and (1,-2). Therefore, we 
distinguish four regions in the plane, as in Figure [SI 



(b) 



P\ 



(a)- 



P*- -tP 



^\rA 


i?4 






. ^2 . \ 


\ 




i?l 






V 


\ •i?^ = 


i?l 




\ 


\^\ 






i?4 


V' \ 



p 

Figure 8. The 5- 
point and 7-point sten- 

cils of system ^M, piGURE 9. Distin- 

and their convex hulls • u j • t 

guished regions tor sys- 
tem (|XTO)) 

For each of the distinguished regions i?' we will project each stencil of the system 
onto a line with direction s g i?*. On these lines we will indicate (schematically) 
the ranges of n- values for both field variables a,p (our ranges start at 0). This 
procedure is similar to the scalar case, where the difference d* is equivalent to a 
range for n. For every equation in a system we need a range (and hence a difference) 
for every field, and also their relative positions. Moreover, once we have this data, 
there is still some freedom left how to position the equations with respect to each 
other. In fact, as we will show in Appendix [Cl we have to split region i?^ into 
two regions where the two equations are positioned differently with respect to each 
other. Also we will see that for regions R^ and R^ we may take the same initial 
values. However, updating the initial values is done differently in those regions. 

For s G i?^ the ranges are given in Figure [TUl 



_aj> 



2a 

—t — 



2a + b 4a + 26 4a + 36 

1 1 1 



Figure 10. s- Reduction for system p.l9|) . with s e N^ 



We now have to find a range of n for each variable a, p such that each equation 
can be used to update one of the variables. Moreover, we want to minimise the 
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number of field values that have to be calculated to update both fields. Therefore, 
we move the first equation in Figure fTOl over a distance b with respect to the other, 
see Figure [TT] 



•- 



-A V 



^. e. fi. 2. 

b 2a + b 2a + 2b Aa + 2b Aa + Sb 

I 1 1 1 1 1 

Figure 1 1 . Alternative s-Reduction for system (|3.19p , with s e N^ 

For system p.l9[) with s e i?^ one can see that if we take the range < n < 
4a + 36 for a and the range 2a + b < n < 2a + 2b for p, then using the equation 
()3.19ap we can determine the values oi p a.t 2a + 2b < n < Aa + 2b. Subsequently, 
the equation (j3.19bl) can be used to calculate the values of cr at ri = 4a + 3b. We 
have obtained the first item in the following proposition, whose complete proof is 
given in Appendix [Cl 

Proposition 3. Let s £ 1? be such that (a, 6) is not equal to (1,0) or (1, —2), and 
let W be as in Figure\^ Then system h3.19]) admits a well-posed s-periodic initial 
value problem. 

• For s e R^ the set 

{crP,pP^ : n E Nia+sb, m - 2a - 6 £ Nb,p £ N J 

provides a well-posed initial value problem of dimension 4| si + S2 |- 

• For s E R^ U R^ (including the boundary between R^ and R^ ) the set 

{aP,pP^ : n e N4a-b,m ENt,pE NJ 

provides a well-posed initial value problem of dimension 4| si |. 

• For s e i?^, b < 3a, the set 

{aP,pP^ : ne N3b-4a,m + 5a - 2b E 'Nb,p E Nr} 

provides a well-posed initial value problem of dimension 4| si + S2 |- 

• For s E R^, b>3a, the set 

K, pP, : n e N3b-4a, m + 4:a-bENb,pE Nr} 
provides a well-posed initial value problem of dimension 4| si + S2 |. 

4.2.1. (2,-l)-reduction of system i3.19fl . We specify initial values between two par- 
allel lines that squeeze the cr-stencil, as in Figure[T2l The initial values are denoted 
crJJ = a:n+i and p^ = Vn+i- Clearly, from these initial values we can calculate the 
values ?/4 and y2, j/i, using equation (|3.19bp . Then we know enough p- values to 
be able to calculate, using equation p.l9ap . xg, which (to our surprise) does not 
depend on j/3. Therefore, the cr-field is determined by a 7-dimensional mapping, 
namely 

Xi 1-^ Xj+i, 1 < 4 < 7, 



X2X-J X2X^Xq X^XiXj X^Xs 

X7 i-s- 



2^2 (4.6) 



xi X1X4 X1X5 xia;4a;5 
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,a;5, ,xe^ ,a;7--«f8 o 



,X3 '^i '^5 ^e jcr 



2^3 

^1 J^l JES. X4 ,X5 
,2:1 . '.■,X2 X3 

m 

Figure 12. (2,-l)-periodic initial values 



We know enough values to determine the monodromy matrix L = Lq{Mq) ^Lq. 
Its characteristic polynomial is given by 

Det(^I -£)= M^((/i - A)2 + (p- A)/i + 1) - A(/iV,/' + fih + 1), 

where 

X3X4 X3XG X^Xtj X1X4 X5X2 XiXj X5X2 



h 



xiXq X4,X5 XTX2 X2X3 X3X4 x^xe xixe 
xexa xlx2 xqxI 



XtX2 XTXiXi XTXiXi 

xiXr — X5X3 xl x'i 



J 

X2XQ X2X4 X^Xq 

The integrals /i, JJ', J + J' are functionally independent. 

4.2.2. (2,-3)-reduction of system Ii3.19\) . We specify initial values between two par- 
allel lines that squeeze the cr-stencil, as in Figure [T51 From these initial values we 



# 


■a, 








^l 


•^4\o 








o\ 


«jf 2 .^^_ 











°\,#. 


\ 








\^1- 


.,a;4\ 


■ 



Figure 13. (2,-3)-periodic initial values 

can calculate the cr- value (xq) on the dashed line. And once we have done that we 
can calculate j/4. We find an 8-dimensional mapping, namely 

Xi >-^ Xi+i, 1 < i < 5, 

, X4a:3 + yiX4 - X3y2 

X5 H- Xr^ = , 

■ 1 o^ (4.7) 

2:43:5 + yix'^ 

2/3 1-^ ^• 

X3 
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One has to calculate quite a few other points to be able to evaluate the mon- 
odromy matrix, whose characteristic polynomial is given bjo 



where 



X3X2 1 ys.yi 1 2/32/2 1 2/22/1 2/1 

XlXi X^X^ X^XA X3X4 x\ 


X5X1 XzX2y2 

xl xixl 


X'iXA X2yi X32/2 ^3X42/3 2/1 

X5X2 XiXa X5X2 X\X2 x\ 


4 


xiXf,y2 ysxi X5X1 X5X3 

x\x2 XiX2 XaX2 x\ 





Det(jul - -£) = /i^((ju - Xf + i-ih + \{JJ' - /i - 1) + 1) - A^, 



/i = - 



J = 



The integrals /i, JJ'^ J + J' are functionally independent. 

4.3. Initial value problems for the QQD scheme. We define regions in the 
plane as follows 

= {seZ^ :a > 0,6>0,e = 1}, 

= {seZ^ :0 < 6< a,e = -1}, 

= {seZ^ :a < 6< 2a,e = -1}, 

= {sgZ^ : 2a < fe< 3a,e== -1}, 

= {seZ^ : 3a < 6,e = -1}, 



R' 
R^ 
i?'5 



see Figure HH We prove Proposition |4] in Appendix [ 
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•i?2 ■ 


i?l 




R^ 
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Figure 14. Distinguished regions for the QQD scheme 



Proposition 4. The QQD scheme IIS. 20\) admits a well-posed s-periodic initial 
value problem if (a, 6) is not equal to (1, 0) or (1, —2). 



■^The factorization oi I2+ Ii + l where I2 is the coefficient of ^^ A is due to Claude Viallet |40| . 
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• For s E R^ the set 

{<, <„, wf -.ne N3a+26, m - 2a e Na+b, I - 2a - b eNb,p eNr} 

provides a well-posed initial value problem of dimension 4| si + S2 |. 

• For s e -R^ the set 

{<, <, wf -.n-be E3a-b, m e Na, leNb,pe N,.} 

provides a well-posed initial value problem of dimension 4| si |. 

• For s e i?'^ the set 

{<, vP„wf -.n + a-be N^a-b, m eNa,le Nb,p G N,.} 

provides a well-posed initial value problem of dimension 4| si |. 

• For s <E R^ the set 

{<,<,u;f : 77, - 3a + 5 G N2b-3a,m e Nb-a,l e Nb,p G NJ 

provides a well-posed initial value problem of dimension 4| si + S2 |- 

• For s E R^ the set 

{<,<,wf : 71 - 2a G n2b-3a,m- a G Nb-aJe Nb,p G NJ 

provides a well-posed initial value problem of dimension 4| Si + S2 |. 

In the following reductions we adopt the notation u'^ ~ a;„ , ?j^ = y„ , ui[^ = z„ . 

4.3.1. (2,-1) -reduction of the QQD scheme. Updating the initial values as described 
in Appendix [D] yields the following mapping 

Xi H> Xi+i, 1 < i < 5, 

, xiyi 

X5 ^ —, 

Vi 
yo^2/i, (4.8) 

, , , x^XiZo 

2/1 ^-> 2/1 = a^i + ?/o - 2:5 - zo H , 

a;ia;2 

Xa^Zq 

Zq ^ . 

Xl 

The characteristic polynomial of the monodromy matrix is given by 

Det(Ml - £) = ^?{{^i - \f + {n- A)/i - XJJ' + h) + X{^xh - h), 
where 

h = (2/0 - zo + Xl) f- yi(zo + X5- yo), 

Xl 

h = z^XiX^yi, 

h = izaX2 - x^yi - X2X1 - yoX2)x3X4 - xiyiix^x^ + X2X5 + X2X3), 

h = xiX2X'iXiXzyi, 

J = Xl +.T3 + J/o - Zq. 

The integrals /i, /2, /a, /4, /s = JJ', Ig = J + J' are functionally independent. 



A HIGHER ANALOGUE OF DISCRETE-TIME TODA AND A QQD SCHEME 23 

4.3.2. (2,- 3) -reduction of the QQD scheme. Updating the initial values as described 
in Appendix [D] yields the following mapping 

Xi^Xi+i, 1 = 1,2, 



X3 I— > x\ 



3 



(4.9) 



yo ^vi 

yi '-^ yo + Z2 - zq, 

Zi^^Zi+i, i = 0, 1, 

The characteristic polynomial of the monodromy matrix is given by 
Det(^I - L) =/i2((^t - A)2 + ^i{{h - /2 + l)/i - /3 - J J') + Kh + (2/2 - 3)/i)) 



where 



yo 






£2 _^ yo _j_ yo ^ y]_ 

a;3 xz xi X2 ' 






X2ZQ - X2yo - X2Z2 + Z1X3 - 


- xzyi - 


-yixi, 


Z0Z1Z2 







+ A^(m(3-/2)-A)/i, 



_ a;2a;ia;3 
^1 — 

h = 

h = 

h = 

X3 

J = yo~ zo 



5. Dimensional reduction 

In this section we present the explicit formulae for the reduced variables, required 
to implement the dimensional reduction of the mappings obtained in section 4 from 
periodic reduction from the HADT equation and from the system p.l9p . These 
formulae arise from the symmetries of the partial difference equations. In certain 
cases an extra scaling symmetry of the integrals has to be used to obtain the 
dimensional reduction. 

Whereas the dimensions of the s-reductions from the HADT equation, of the 
intermediate system, and of the QQD scheme all coincide, these reductions differ 
both with regard to the number of integrals we found in the previous section and 
with regard to the number of symmetries we will encounter in the present section. 
The question whether the s-reductions of the various equations are related, e.g. 
through a change of variables, is left open. 

5.1. The HADT equation. Equation (|3.18p is invariant under the following con- 
tinuous symmetries: u >—^ ue and u i— > ue^^^' . The first yields a symmetry of 
the mappings, whereas the second gives us a 2-symmetry of the mappings. Their 
generators are 

^ Q ^ ^ 

i=l 4=1 



(5.1) 
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A set of joint invariants of these vector fields is given by 

, XiXi X2X5 X3XG XaX-j X^Xs XqXi -, 
{Wi = ,W2 = ,W3 = ,W4 = ,W5 = ,We = }. 

2;2a^3 2;3a;4 x^x^ x^xq xqX7 xrxg 

In these variables the s — (2, —1) mapping (I4.3P reads 
Wi i-> Wi+i, 1 < i < 5, 
/ 

W5 *-^ , f = W2 - W5 + {Ws - W2)W4 + W1W2W3W4W5, 

'W2W3W4'W5 

W2W4'We 

which has three functionally independent invariants 

/i = wi + W5 + u'2W4(l — W1W3W5) H , (5.2a) 

'W2W3W4 

_ {W1W2W3W4 +W2+W3- 1){W2W3W4W5 +W3+W4-I) , , . 

W2W3W4 

I3 = W2 + W3 + W4 — W4W2 + W1W2W3W4W5. (5.2c) 

The 2-integral J is scaled by the generator of the second symmetry and so it does not 
survive the reduction. The s = (2,-3) mapping ()4.4[) reduces to the six dimensional 
mapping 

Wi ^ Wi+i, 1 < i < 5, 

9 1 1 

W5^ , g^WiW2W3'i /t- o\ 

W4W5 W3 W4 W-'J,' 

W4WQ 

We h-^ 

wiwsg 

and has three functionally independent invariants 

11 = W3{wiW2 + W2W4 + ^4^5) H , 

W3 

12 = W2W3W4{wiW2 + WlW^ + ^4^5) — ^2^3 — W3W4 + W2W4, 
J J' = W3{wiW2W3W4 — 1){W2W3W4W5 — 1). 

5.2. The lattice system (j3.19p . For the case s = (2,-1) the mapping has a 
symmetry and a 2-symmetry, which arise from scaling symmetries of the system. 
Their generators are given by 

i—i i—i 

The integrals admit an extra scaling, which is neither a symmetry, nor a /c-symmetry, 
of the mapping, namely 

d d d d ^ d 

axi 0x2 0x5 oxq oxr 

Still, we can use joint invariants of these three vector fields, 

a;ia;4 X5X2 X3Xe X4Xj 

Wi = ,W2 = ,W3 = ,W4 = , 

^2X3 a;3a;4 X5X4 X5Xe 

to reduce the dimension of the mapping by 3 dimensions. This is due to the fact 
that if we reduce the third vector field by the first two the resulting vector field is a 
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symmetry of the reduced mapping. The mapping (|4.6p reduces to a four dimensional 
mapping 

Wi t-^ w,;+i, 1 < i < 4, 

W2W4 — W2 — w^ — Wi (5.4) 

Wi M- , 

WIW2WZWA 
which admits two functionahy independent integrals 

WlW^ + W2W4 — U'l — W2 ~ w^ ~ W4 



/l = Wi + W2 + W3 + W4 + 



JJ' ~ {Wi + W2 + W3 + W4 — W1W3 — W1W4 — W2W4) 



W1W2W3W4 

W1W2W3W4 + W2 + W^ — 1 



WlW2W3Wi 



(5.5) 



Similarly wc find three reductions for the s ~ (2, —3) case. Here choosing vari- 
ables 

yi 2/2 2/3 X2X3 x-iXi 

Wi = ,W2 = ,W3 = ,W4 = ,W5 = 

Xz X4 X5 X1X4 X2X5 

reduces the mapping (j4.7p to the five dimensional mapping 

Wii^Wi+i, i = 1,2,4, 

UI3 l-> U'l + UI5 , 

1 

«75 h^ Ug = -— . 

(f + Wi - W2)W4W5 

This mapping admits the two integrals 

/l = Wi{w2 - Wi) + W2iw3 - W2) + Wsiwi - W3) + W4{W2 - Wi - 1) 

I 

+W5{W3 -W2-I) , 

w^w^ 

J J' — {1 +W2 - W3 - W4)( 1 - W2 + W3 - W5 - (Wi — W3 + W5){wi - W2) 

W4W5/ 

We note that in both cases the function J is no longer a 2-integral of the reduced 
mapping, and that J + J' can not be expression in the reduced variables. Also, 
in the s = (2, —3) case wc have only two integrals for a five dimensional mapping 
which is not enough for complete intcgrability. The reduction was done using 3 
(scaling) symmetries of the integrals, but the function J is only invariant under one 
of them. It turns out it is also invariant under a linear combination of the other 
two. Using this, we can then perform 2-reduction for mapping (j4.7p . retaining all 
three integrals. Choosing variables 

xl x\ x\ y2 ys yixz 

Zl = ,Z2 = ,Z3 = ,Z4 = ,25 = ,Zq = ^, 

xsxi a;4a;2 a;5a;3 X4 X5 2/2 

the mapping reads 

Zi i^ z^+i, i ^ 1,2,4 

23 i-> 4 == zizlzl{l + zqzIzs - Z4) 

Z5^^Z3(4 + ■^4■^6) (^-^^ 

Zi 

ze i-^ ~r3' 
^3^5 
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which has three invariants: 

Z6Z3zl{zezlz3+ziZ2-Z4-Z5)+zl+z'^-Z4Z5-\ 2 (24-1)212:2 -(25-24-1)2:322, 

212223 

and J + J' and JJ', where J = (1 + 24 — 25 — 2i22)(2i2223)~^. 
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Appendix A. Some determinantal identities 

In the establishment of the recursive structure we need a number of determinantal 
identities, which we derive using the Sylvester Identity. So we present a proof of 
the Sylvester identity, which was first presented by Kowalewski [20] , Bareiss [2] and 
Malaschonok [241 ESj and these seven proofs are presented together in [1] . 

We consider an (n + m) x (n + m) matrix R with elements r^ and determinant 
i?|, also written Det(i?). Then we partition R and factor by block triangularization 
such that 



R = 



A B 
C D 



A 
C 1 



1 A-^B 

D-CA-^B 



where A is a nonsingular square matrix of order n, then 

\R\ = \A\.\D-CA-^B\. 
If we multiply both sides by |^|"'^^, this becomes 

\Ar-^\R\ = \\A\{D-CA-^B)\ 



(A.l) 



(A.2) 



because the determinant on the right side of (jA.2[) is of order to. We can reduce 
this equation further to 



1^ 



rn— 1 



i?| = \\A\D-CABl 



(A.3) 



since A ^ ~ t4t (where A represents the adjugate matrix of the inverse matrix 
A^^), and the determinant of A is assumed to be ^ 0. Specifying some entries in 
(|A.1I) . taking ^ to be an rt x n block and D to be an ?7i x jti block, we have the 
formula: 



(Det(A)) 



7n—l 



A 



"-1 



bi 



D 



Det, 



iXTTl S 



-L^^tn X n (-^J -^f 



{clAbA 



in which the full matrix is supplemented with to n-component column vectors hi 
and m n-component row- vectors c*. If we consider the case tti = 2 ie. the removal 
of two rows and columns, then we get then determinant identity 



(HetiA)) 



- + 
4 



bi 


&2 


dii 


di2 


d2i 


1^22 



= Det2 



(Det(A)) 



dii di2 

d21 1^22 



c{Ah-i c{Ab2 
C2Abi C2Ab2 



[Det{A)dii ~ clAbi] [Det(A)d22 - C2A62] 
-[Det(yl)d2i - ciAb-i] [Det(A)di2 - 0*1^62] 
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which can be symboUcaUy written as: 



(A.5) 



(where the red hnes denote rows and columns omitted from the original determi- 
nant). It is then necessary to reorder the position of the row and column to tailor 
the identity to our requirements. 

While (|A.5I) is the key identity by which the recurrence structure for ordinary 
one- variable orthogonal polynomials is obtained, for the elliptic two- variable orthog- 
onal polynomials we need (in addition to (|A.5|) ). dcterminantal identities involving 
the simultaneous removal of more than two rows and columns. Thus, the main 
identity we use from the general formula (jA.4|) will be the case to = 3, leading 
to the different recurrence relations for p.3|) and p.7p . So considering m = 3 we 
obtain from (|A.4p the following 3-row/column Sylvester type identity: 























X 






















(A.6) 

which is obtained from the expansion of (|A.4p for ttt, = 3 after recombination of 
terms using the earlier 2-row/column Sylvester identity (|A.5[) . 
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Appendix B. Intermediate determinants 

Previously wc have chosen to apply a 3 row/column Sylvester identity to the Pf. 

and QJj, polynomials since this identity (JA.6P does not introduce new determinants. 
However for some 2 row/column Sylvester identities it is possible to control new 
determinants, which can be removed to still give equations in term PJ: and Qj, 
only. These " intermediate" determinants are introduced here. 

B.l. The linear Pf, equation. By considering the PJ: with the first column and 
penultimate row removed we get the intermediate quantity: 



Tl!lA^,y) 



(ei,e2) 
(e(+i,e2) 



(ei,e3) 
(e;+i,e3) 



(e;+fc-3,e2) (e;+fc-3,e3) 

62 63 



together with a corresponding Hankel determinant: 



(6;,efc) 
(6; + i,efc) 



(6(+fc-3,6fc) 

efe 



/ni?l2 



(B.l) 



n 



(0 



(6(,62) 
(e;+i,e2) 



(6;, 63) 

(6i+l,63) 



(ei,6fe) 
(e;+i,efc) 



(B.2) 



(ei+fe_2,e2) (ei+fc_2,e3} [ei+k-2,ek) 

Using the usual Sylvester identity we can now derive the following two equations 





n!" 




p(0 







p, 



(0 



T, 



C+i) 



fe-i 



A('+i)n(') 

^fc-2 '^^k-l p(l+l) 
^fc-l-^^fc-2 



a(0 tt('' 

p{l)_rj.{l) ^fc-2%-1 p(0 

^fc-l^^fc-2 



(B.3a) 



(B.3b) 



Eliminating the T^ polynomials in favor of the PJ^ polynomials and using a Hankel 
identity 




hWaC+s) 



n(0 A ('+3) 



(0 






(B.4) 



allows the derivation of the linear relation in Pt p.l2p 
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B.2. The bilinear relation in Aj^ ,Qj. . It is possible to combine some three- 
term bilinear Hankel relations, to give a three-term relation in terms of Aj. , 8^ . 
This particular identity is achieved with the introduction of a new intermediate 
determinant S which is essentially Q with the first column and last row removed. 



A^U^l? =. Q^Lu2ll' aL+:'s('L (B.5a) 





A« 


et" 



'^k-2 



■'k-l^'-k-l 



{l)jj{l + 2) _ cv(0 yr{l + 2) 









'■k-l "fe-1 



,('+2)^(0 



(B.5b) 



where 



^(0 



(ei,e2) (e(,e3) 

(ei+2,e2) (e;+2,e3} 



(e(,efe) 
{ei+2,ek) 



(B.6) 



S is then eliminated in the combination of the two bilinear relations (jB.Sp . 
This equation can be expanded (in the third term) using a third bilinear equation 




^fc ^^fe-2 ^ ^fe-l^^fc-l ^fe-l ^^fe-l • 



(B.8) 



Thus we are left with an equation in terms of A and Q only. 



,('+!) /o(0 rr('+2) 



Ar;^(e^inri^-ernr 



(l)jj{l+2)s 



(l+l)jj{l+2) 



e^,(Ar^^nrr^-Ar^^nno 



(^+l)TT('+2)^ 



A('+2)A(0n('+2) 

"^fe-l ^fc ^^fc-2 



After we have canceled the necessary terms, we are left with (|3.17bl) . 



B.3. Relations in P^ and Qj, . For relations involving the Hankel determinants 
Pf. and Qj, we can derive a recurrence relation xQj. and two linear relations. By 
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implementing a cutting of rows and columns on the matrix for P^, according to: 



and applying the 3-row/column Sylvester identity in that situation, we obtain: 



Pt'' - -(t^ 



a(0 C)('+2) 
'-^k-2^k-2 p{l) 

aW o('+2) fe-1 
^fe-l^/c-3 



C)(0 a('+2) 



Qi'li , ?^0,l. (B.9) 



For the linear relations we first obtain 




A 



('+2)p(0 



-,(0 _ rp{l+2) _ ^fc-2 ^ fc-1 p(i+2) 
^fc-l^^fc-2 



(B.IO) 



which by eliminating the T^ ' using (|B.3aP together with the Hankel identity (jB.5a[) 
leads to 



(0a('+2) 



(^=Pl 



A.'^A 









(B.ll) 



i« 



3(0 ;c 



This three term equation in terms of Q), and P^ is similar to p.l2p . To find an- 
other of this type of equation we must first introduce the intermediate polynomials 

oil) . 



stU^^v) 



{ei,e2) (£(,63) 

(e;+2,e2) (e;+2,e3) 



{ei+k~2,e2) (e;+fe_2,e3) 

62 63 



(ez+2,efc) 



(ei+fc-2,efc) 



/ri?l2 



(B.12) 
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Then applying the followmg Sylvester identities to PJ: and Q\. , we subsequently 
obtain the relations: 



Qr 



(0 



,(0 



^^' 



g{l) ^ _ ^k-l'^k-2 Q(l) 



'- fe-2^fc-l 



P, 



(0 



q{l) 



^^k-l^k-2 ^(l) 
-L k-2'-^k~l 



We eliminate S'j,_j^ and making use of the Hankel determinant identity 



^ ^fc "fc-2 ~ '-^k-l^k-l ^k-l^^k-l 




which leaves: 



Qi'^ 



p: 



(i) 



(Oo(0 



A:'e 






fc-2 ^(0 



(B.13a) 



(B.13b) 



(B.14) 



(B.15) 



A HIGHER ANALOGUE OF DISCRETE-TIME TODA AND A QQD SCHEME 33 

Appendix C. s-Periodic initial value problems for system (|3.19p 

Proposition |3] is obtained by considering the ranges of the n- values for the equa- 
tions with periods taken in the remaining regions. When s G i?^ (see Figure IH]) we 
get the following reduction. 

P. ^^ ^ 



a^ 



-^ 



h 2a ~b 2a 4a - 6 

I 1 1 1 1 

Figure 15. s- Reduction for system (|3.19p . with s e i?^ 

So we first need to calculate the values oi p ior b < n < 2a — b. This can be done 
using equation (|3.17ap because the values of cr at ri smaller than 2a + {2a — b — b) = 
4a — 2b are given initially. Then the values of cr at n = 4a — 6 are determined by 
(j3.17bp . When s £ i?^ we get 

P. ^ fi. ^ 



2a — b a b a + b Aa~ b 



Figure 16. s- Reduction for the system p.lQp . with s e i?^ 

In Figure [Tni the dashed line indicates that a + b can either be to the left, or to 
the right of 4a — b. We have a -I- & < 4a — 6 in the region 3a > 26 > 2a. Here it does 
not matter whether one first calculates the values of p bX n = b, using equation 
p.l7ap . or the values of tj at n = 4a — 6, using equation (j3.17bp . In the region 
3a < 26 < 4a we have a + 6 > 4a — 6. Here we need first calculate the values of a 
for 4a — 6 < n < a + 6, using equation (j3.17bp . This can be done indeed, because 
the values of p with n < 2a — b + (a + b) — (4a — 6) ^ b — a are given initially. 
Subsequently the values oi p at n = b can be obtained. 

When s G i?^ we also need to distinguish two cases, but here this leads to 
different initial value problems. When 2a > 6 < 3a we can first calculate the values 
of (T at n = 36 — 4a, and then the values of p at n = Sb — 5a, see Figure [TT] 

p a p a 



\ e, e. ^ 

26 -5a 6 -2a 26 - 4a 36 - 5a 36 - 4a 



Figure 17. s- Reduction for the system p.l9p . with s G i?*^, 6 < 3a 

When 6 > 3a we need first calculate the values of p at n = 26 — 4a, and then the 
values of (7 at n = 36 — 4a, see Figure [HI 
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(J n p a 

6 -4a 6 -3a 6 - 2a 26 - 4a 26 - 3a 

f 1 1 1 1 1 1 

Figure 18. s- Reduction for the system p.l9p . with s £ i?*^, 6 > 3a 



We leave it to the reader to check that the given initial value problems can also be 
updated in the negative n-direction. In Appendix |E] we illustrate with an example 
that this is not necessarily the case, and thus has to be verified separately. 
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Appendix D. s-Periodic initial value problems for the QQD scheme 
The stencils of the QQD scheme (|3.20|) arc depicted in Figure fT9l 





w 



u- w 



Figure 19. Graphical representation of the QQD scheme 

To describe initial value problems for the QQD scheme we follow the following 
procedure. We project the stencils of the three equations onto lines with directions 
in different regions and translate them in such a way that, given the values of the 
fields for certain ranges of n, each equation determines a value for an update of one 
of the fields. The proof of Proposition consists of 5 pictures, one for each region. 
For every picture one has to check the following: 

(1) that the ranges of n-values for u, v, w correspond to each equation in the 
system, for any chosen direction s in the particular region, this includes 
checking the order of the linear expressions in a, h. 

(2) that, given the initial values for u,v,w, every equation can be used to 
update one of the fields (in both directions) . 



*, 






v_w 




v_w 


"■ 


i_ 




V 






U V 






"• 




w 


UJU 








1 


b 

1 


2a 

) 


2a + b 

1 


2a + 2b 
1 


3a + 6 

— J 


3a - 

1 



Figure 20. s- Reduction of the QQD scheme, s G i?^ 



We will perform check (1) and check (2, positive direction) for the first picture, 
and check (2, positive direction) for the other pictures. This enables one to perform 
s-reduction for any given s G UiW. The rest of the proof is left to the reader. 

When e = 1, as in region R^, with one step to the right (downwards) n increases 
by b (by a). Let n be equal to on a line with direction s G R^ through the first 
point u in the first equation of the system (the upper-left point in Figure I19p . We 
let the upper-left point from the second equation coincide with the one of the first 
equation and we let the upper-left point from the third equation be one step to the 
right, so that it is at distance b from the others, as in Figure [^ Moving the line 
over the Figure [3T] in downward direction, depending on the direction s, the line 
will either first cross this point at distance b, or the point in the second equation 
at distance 2a. This is indicated in Figure [20] by the dashed line between b and 2a. 
Next the line moves over the point with distance 2a -I- b, and then it depends on s 
whether 2a + 26 or 3a + 6 is encountered first. This is again indicated by a dashed 
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Figure 21. The equations of the QQD scheme on top of each 
other, and a hne with direction (2, 1) G R^. 



line. Thus, Figure [20] represents the ranges of the n- values of the three fields of the 
three equations in the QQD scheme and their relative positions, for s S i?^. 

Next, suppose initially we are given the values of m at < rt < 3a + 26, of w at 
2a < n < 3a + b, and of w at 2a + 6 < n < 2a + 2&. Then, for all s G 7?^, we can 
use the second equation to determine the value of w at n = 3a + 6. If a < &, then 
we know all values of w at 2a + 6 < n < 3a + 6, and hence we can use the first 
equation to calculate the value of m at n = 3a + 26. Then, we can calculate the 
value of w at n = 2a + 26. However, if a > 6, then we first have to use the last 
equation to calculate all values of w at 2a + 26 < n < 3a + 6, before we can use the 
first equation to calculate the value of u at ti = 3a + 26. In both cases, the initial 
value problem is well-posed. 



u V 



uw 
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3a — b 



3a 



Figure 22. s-Reduction of the QQD scheme, s e _R^ 

In R^ we first use the third equation to determine the values of w at 6 < n < a, 
then the first equation to determine the values of w at n = a, and finally the second 
equation to determine the values of u at n = 3a. 
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Figure 23. s-Reduction of the QQD scheme, s e R^ 
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In R^ the values of u at n = 2a and the values of w at n = a can be determined 
independently, using the second and first equation, respectively. After having de- 
termined the values of w at n = 2a we are able to find the w at n = b, using the 
third equation. 

V W U V w u 



V U V u 



u uw 



3a — 5 a b — a 2a b 

I 1 1 1 1 1 

Figure 24. s-Reduction of the QQD scheme, s £ i?"^ 

In i?^ the values oi u at n = b and the values oi v at n = b~a can be determined 
independently, using the first and second equation, respectively. After having de- 
termined the values oi u at n = b we are able to find the w at n = b, using the 
third equation. 
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Figure 25. s-Reduction of the QQD scheme, s G i?^ 

In R'^ the values of v and w at n = b can be determined independently, using 
the second and third equation, respectively. Next the values of u at n = 25 — a can 
be determined using the first equation. 



38 PAUL E. SPICER, FRANK W. NIJHOFF, AND PETER H. VAN DER KAMP 

Appendix E. A non-invertible reduction of the QQD scheme 

Using the same notation as in scction l473j consider initial values {xq, xi, X2, 2:3, X4, 
yo,yi,zo, Zl} for the (2, —1) reduction of the QQD scheme. We use equation p.20b[) 
to determine X5 , equation (j3.20cp to determine Z2 and equation p.20a|) to determine 
2/2- Thus we find the nine dimensional mapping 

Xi 1— > Xi+i, < i < 5, 

yoxo 



(D.l) 

Zq + Z2 — X5, 



X4 


K^ 


X5 




yo 


ya 


n> 


yi, 






2/1 


1— >■ Xi 


+ 


yo ~ z 


zo 


^ 


Zl, 






Zl 


h^ 


Z2 


= 


X5 

Zl. 

X2 



for which the staircase method docs yield five functionally independent integrals. 
However, the mapping is not invertible and so the initial value problem is not well- 
posed. This example illustrates that in order to prove the well-posedness of an 
initial value problem one has to show the initial values can be updated in both 
directions. 
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